NASA Technical Paper 1367 


NASA 

TP 

1367 

c.l 



An Analytical Technique for ; 

Predicting the Characteristics 

of a Flexible Wing Equipped With 

an Active Flutter-Suppression System 

and Comparison With Wind-Tunnel Data 


Irving Abel 


FEBRUARY 1979 





J 


TECH LIBRARY KAFB, NM 




NASA Technical Paper 1367 


tech library kafb nm 

l mill HI.. r °» 



013437a 


An Analytical Technique for 

Predicting the Characteristics 

of a Flexible Wing Equipped With 

an Active Flutter-Suppression System 

and Comparison With Wind-Tunnel Data 


Irving Abel 

Langley Research Center 
Hampton, Virginia 


NASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
information Office 


1979 


SUMMARY 


An analytical technique for predicting the performance of an active 
flutter-suppression system is presented. This technique is based on the use 
of an interpolating function to approximate the unsteady aerodynamics. The 
resulting equations are formulated in terms of linear, ordinary differential 
equations with constant coefficients. This technique is then applied to an 
aeroelastic model wing equipped with an active flutter-suppression system. 
Comparisons between wind-tunnel data and analysis are presented for the wing 
both with and without active flutter suppression. Results indicate that the 
wing flutter characteristics without flutter suppression can be predicted quite 
well but that a more adequate model of wind-tunnel turbulence is required when 
the active flutter-suppression system is used. 


INTRODUCTION 

A difficulty in analyzing active flutter-suppression systems lies in the 
modeling of the unsteady aerodynamic forces. These aerodynamic forces are nor- 
mally computed only for simple harmonic motion at discrete values of reduced 
frequency. The use of harmonic motion for flutter analysis is adequate since 
the problem is one of finding the neutral stability boundary for which the 
motion continues with constant amplitude. The problem facing the analyst is 
one of modeling the unsteady aerodynamics for arbitrary motion. 

In lieu of developing a completely new aerodynamic theory, there has been 
considerable interest in using the results of oscillatory unsteady aerodynamics 
to generate approximate solutions for arbitrary motion (refs. 1 to 4) . This 
paper presents a method for analyzing active flutter-suppression systems. The 
method is based on a technique for approximating the unsteady aerodynamics in 
the time plane through an interpolating function in the frequency plane. By 
using the aerodynamic approximating function, the equations of motion are for- 
mulated in terms of linear, ordinary differential equations with constant coef- 
ficients. Active control functions are added to the equations in a straight- 
forward and convenient manner. The resulting equations are reduced to a series 
of first-order differential equations which are solved to construct a root 
locus of the modes as a function of dynamic pressure. Also included is a 
method for calculating the response of the control system to turbulence. 

The analytical method is then applied to an aeroelastic model equipped 
with an active flutter-suppression system. Comparisons between wind-tunnel 
data and analysis are presented for the wing flutter characteristics both 
with and without active flutter suppression. Also presented is a comparison 
of wind-tunnel data and analysis for the response of the active flutter- 
suppression system to tunnel turbulence. 



SYMBOLS 


Ai(s) polynomial in s (see eqs. (B3) ) 
a speed of sound 

b reference semichord used in aerodynamic theory 

c streamwise local chord 

D(s) denominator polynomial in s of transfer function 

1^ feedback filter parameter 

E error function 

H ( 03 ) control surface frequency response function 

j -fT 


k 

L 

L s 

1 

M 


Mi 

m(x,y) 

N(s) 

n 

Ap(x,y,t) 


Qi 


Qij 


reduced frequency, wb/V 

characteristic length in Von Karman gust spectrum 
wing span 

Laplace transform operator 
Mach number 

generalized mass in ith vibration mode 
mass distribution 

numerator polynomial in s of transfer function 

number of flexible modes 

pressure distribution 

generalized aerodynamic force in ith mode 

generalized aerodynamic force in ith structural mode due 
to pressure distribution in jth mode 


q 


dynamic pressure. 


-pv 2 
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q^ generalized displacement in ith mode 

r number of control surfaces 
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S integration surface 

s = jw 


t 

V 


w g 

x,y 

Z(x,y) 

h 

*a 

«c«c' 

c 

p 



4>i(io) 

4>i (x,y) 
a> 


“n,i 

Matrices: 


time 

free-stream velocity 
gust velocity 

streamwise and spanwise coordinates, respectively 
vertical deflection 
aerodynamic lag 
control-surface deflection 

control-surface command and compensation, respectively 

viscous damping coefficient 

fluid density 

rms gust velocity 

Von Karman gust spectrum 

normalized modal deflection in ith mode 

circular frequency 

circular frequency of ith natural mode 


[A] 

[Ac! 

[Ail 

[Ai,6l 

{A if g} 

[a]. . 
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matrix representing first-order equations of motion, system off 
matrix representing first-order equations of motion, system on 
real aerodynamic matrix coefficients 

real aerodynamic matrix coefficients for control surfaces 
real aerodynamic matrix coefficients for gust forces 
Ag] see equations (B3) 

approximate aerodynamic gust vector 
real coefficients of equations of motions 



[K] generalized stiffness matrix 

[M] generalized mass matrix 

[N] numerator polynomial in s of transfer function matrix 

[ Q] matrix representing generalized aerodynamic forces 

A 

[Q] matrix representing approximate aerodynamic forces in the 

Laplace plane 

{<?}/{ <j) complex response vectors 
[T] transfer function matrix 

{X} response vector of first-order equations of motion 

[<}) t ] matrix of modal deflection at sensor location 

Subscripts: 
c control 

I imaginary part of complex value 

max maximum 

R real part of complex value 

rms root-mean-square value 

t sensor location 

Dots over symbols denote derivatives with respect to time. 


DESCRIPTION OF TECHNIQUE 

Three-dimensional unsteady aerodynamics are normally computed, at a given 
Mach number, for simple harmonic motion at specific values of reduced frequency 
k. The control law for an active flutter-suppression system is usually given 
as a transfer function which relates control-surface motion to wing response. 

It is normally expressed as a ratio of polynomials in the transform variable 
s. The problem associated with the analysis of an active flutter-suppression 
system is developing a set of equations where the form of the unsteady aerody- 
namics and the control law are compatible. The approach taken in this report 
is to permit the variation of the aerodynamic forces with frequency to be 
approximated by a rational polynomial in the variable s. This technique is 
similar to that described in references 3 and 4. 

The generalized aerodynamic forces are approximated in the s-plane through 
an interpolating function of the form 
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6_ [Ajs 


( 1 ) 


[Q] = 


[A 0 ] + [A T ]( % + tA 2 ](-) s 2 + 

\ v / m=3 


V 1 

3 + i s ”- 2 y 


where s = jo) and 3 m _2 = a) m-2^/ v * As described in reference 3, the form of 
equation (1) permits an approximation of the time delays inherent in unsteady 
aerodynamics subject to the following requirements: complex conjugate sym- 

metry, denominator roots in the left-hand plane, and a good approximation of 
the complex aerodynamic terms at s = jo). The approximating coefficients 
(Ag, A], . . . , Ag) in equation (1) are evaluated by a least-squares curve 
fit (described in appendix A) through the values of complex aerodynamic terms 
at discrete values of reduced frequency. Figure 1 illustrates a typical fit 
through the values of ccxnplex aerodynamic coefficients. The solid curve 
represents the approximating function. 


Wing Without Flutter-Suppression System 

The equations of motion are formulated through a modal approach using 
Lagrange's equations of motion. In the modal approach, the elastic defor- 
mation at any point on the wing is described by a linear combination of 
orthogonal modes. 


n 

z(x,y,t) =!>Iqi(t) <l>i(x,y) 
i=1 

where <|>i(x,y) are the undamped natural modes of the system and n is the 
number of modes used. Assuming a viscous form for structural damping, the 
equations of motion become 


Mi qi(t) + 2? 1 M i a> Ilfi 4i(t) + Wn,i M i = --PV 2 Qi<t) 


(i— 1, 2, . . ., n) 


( 2 ) 


where 


If 


Mi = \ \ m(x,y) 4>i 2 (x,y) dx dy 


is the generalized mass and 


0l <t, - ^ 

S 


Ap(x,y,t) 

q 


4>i(x,y) dx dy 


5 



is the generalized aerodynamic force. The total pressure distribution 
Ap(x,y,t) can be expressed as the sum of the contributions due to each 
flexible mode. Therefore, 

n 

Ap(x,y,t) = Apj (x,y) q-j(t) 

j=l 


where Apj(x,y) is the lifting pressure at point (x,y) due to wing motion 
in the jtn flexible mode. Substituting this expression for Ap(x,y,t) into 
equation (2) results in 


Mi gi (t) + 2CiM iWnfi q ± (t) +a) n>i Mi q ± (t) 


1 n 

2 j=1 


^ (t) JJ 


Cf' APj( x »y) 


<i>i(x,y) dx dy 


(3) 


where i = 1, 2, . . ., n. By taking the Laplace transform of equation (3), the 
equations of motion can be written as 



[ 2940 ^]! 


+ -Wz([Q] ) 

2 


+ t K] { q} = 0 


(4) 


where I is the Laplace transform operator, and 





Apj (x,y) 

q 


s 


<!>i (x,y) 


dx dy 


k 


ij 



U = j) 
(i t j) 


By substituting equation (1) into equation (4) and following the procedures 
outlined in appendix B, the equations of motion are reduced to a series of 6n 
first-order equations of the form 

s{ X} = [ A] {x} (5) 


The eigenvalues of [A] are the roots of the characteristic flutter equations. 


6 



Wing With Flutter-Suppression System 

The equations of motion are formulated in the same manner as described 
for the case with no control. For the wing equipped with an active flutter- 
suppression system, the equation which corresponds to equation (4) is 



where 

[Mel 

[K c ] 

[Q c l 


control- surf ace inertial coupling 
control-surface stiffness coupling 



Ap,q c (x,y) 

q 


4>i (x,y) 


dx dy 


( 6 ) 


For a single control surface and a single sensor (accelerometer) the control- 
law transfer function can be assumed to be of the form 


*( 6 a> N(s) 

2 (z t ) °( s > 

where 

2 (z t ) = s 2 l4> t J {q} 

and 

6 a = q c 


Therefore, equation (7) can be written as 


Sc = 


s 2 n (s) L J 
D(s) 


{q} 


(7) 


( 8 ) 


where 

N(s) numerator polynomial in s 

D(s) denominator polynomial in s 


and I4> t ] is a row matrix of modal deflections at the sensor location. Substi- 
tuting equations (8) and (1) into equation (6) and following the procedures 
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outlined in appendix B, the equations of motion are written as a series of mn 
first-order equations of the form 


s(x} = [AqHx) 


(9) 


where 


m 

6 + highest order term in 

n 

number of modes 


The eigenvalues of [a^ are the roots of the characteristic equation for the 
wing equipped with an active flutter-suppression system. 


APPLICATION OF TECHNIQUE 

To evaluate the adequacy of the analytical method in predicting the per- 
formance of an active flutter-suppression system, stability calculations were 
made for an aeroelastic wind-tunnel model equipped with an active flutter- 
suppression system. Tests were performed in the Langley transonic dynamics 
tunnel. A photograph of the model mounted in the wind tunnel is shown in fig- 
ure 2. Model geometry is given in figure 3. 

As input to the analysis it is necessary to determine a set of generalized 
masses, mode shapes, and natural frequencies of the model. The first 10 elas- 
tic modes were used for analysis purposes. The modes were determined using a 
finite-element model of the wing. The modes cover a frequency range from 
5.23 Hz to 118.15 Hz. Generalized masses and frequencies are presented in 
table I. 


Aerodynamic Properties 

The aerodynamic terms appearing in equation (2) were calculated using 
doublet-lattice aerodynamics by a numerical method similar to that described 
in reference 5. To calculate the pressure distribution on an oscillating wing 
undergoing simple harmonic motion, the lifting surface is subdivided into an 
array of trapezoidal boxes arranged in strips parallel to the airstream as 
shown in figure 4. The lifting surface is then represented by a lattice of 
doublets located at the quarter-chord of each box. The downwash boundary con- 
dition is satisfied at the three-quarter-chord of each box. The downwash is 
computed from the slope and deflection of each structural mode. The lifting 
surface was divided into 210 boxes arranged in 30 streamwise strips with 7 boxes 
per strip. Oscillatory aerodynamic forces were calculated at eight reduced fre- 
quencies (k = 0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.3, and 1.8). 

Each of the aerodynamic terms was approximated in the s-plane through the 
use of equation (1). The B m _2 terms were arbitrarily selected to be 0.2, 0.4, 
0.6, and 0.8, respectively. Figure 5 shows a comparison between the oscillatory 
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data and the approximating function for Q] ■) and Q 22 at M = °- 90 * These two 
aerodynamic terms were selected because of their importance in the flutter cal- 
culations. As indicated in figure 5, the fit is good. In general, errors in 
the aerodynamic forces were less than 10 percent and in most cases less than 
1 percent. 


Control Law 

The active flutter-suppression system that was implemented on the model is 
illustrated in the following sketch: 



where 

6 C ' (s) 4. 37 x 1 0 11 s (s 2 + t^s + 1.169 x 10 4 )(s 2 + 22.7s + 1 432) (s 2 + 21 6.2s + 4.675 x TO 4 ) 

s 2 z t (s) (s + 1 .1 ) (s + D),,) (s + 21 ,6) 2 (s 2 + 43.3s + 2922) (s + 432.4) (s + 486.5) 2 (s + 628. 3) 2 

in degrees per g unit, 

6 c< s ) 2. 795 (s 2 + 179.4s + 8.945 x io 4 ) 

fi c' (s 2 + 350s + 500 2 ) 


in degrees per degree. 


6 a( s > 3.057 x 10 13 

5 c( s ) (s + 214) (s 2 + 179.4s + 8.945 x 10 4 )(s 2 + 747.9s + 1.597 x 10 6 ) 


in degrees per degree, and 

z t = z(x = 0.60c, y = 0.92L S ) 


Therefore, 


5 a ( S ) 3.734 x 1 0 25 s 3 (s 2 + Cfos + 1 .169 x 10 4 )(s 2 + 22.7s + 1 432) (s 2 + 216.2s + 4.675 x 1 0 4 ) 

z t( s ) (s + 1 .1) (s + D)„) (s + 21 . 6) 2 Cs 2 + 43.3s + 2922) (s + 432.4) (s + 486.5) 2 (s + 628.3) 2 (s 2 + 350s + 500 2 ) (s + 21 4) (s 2 + 747.9s + 1 .597 * 1 0 6 ) 


( 10 ) 
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in degrees per g unit. The control surface has a 20-percent chord and is 
located between span stations y = 0.763L S and y = 0.893L S . Locations of 
the control surface and of the feedback accelerometer are shown in figure 3. 

The feedback filter parameter Djp varies with Mach number M and dynamic 
pressure q in the following manner: 

Djn = -83 . 54q - 900M + 1540 (11) 

RESULTS AND DISCUSSION 
Wing Without Flutter-Suppression System 

The eigenvalues of equation (5) are the roots of the characteristic flut- 
ter equation for the wing without flutter suppression (system off) . Since the 
matrix [a] varies with dynamic pressure, a root locus illustrating the varia- 
tion of the flexible mode eigenvalues with dynamic pressure can be constructed 
for each Mach number. A typical root locus at M = 0.90 is given in figure 6. 
(It should be noted that extra roots associated with the aerodynamic poles in 
eq. (1) are calculated when solving for the eigenvalues of eq. (5), but are not 
presented in fig. 6.) Root loci for each flexible mode are indicated in the 
figure. Arrows indicate increasing dynamic pressure. A classical flutter 
behavior is apparent since the frequencies of modes 1 and 2 tend to coalesce 
with increasing dynamic pressure as mode 1 crosses into the unstable region. 

The value of dynamic pressure q at flutter is given in the figure. Calcula- 
tions performed at M = 0.60, 0.70, and 0.80 show a similar behavior but with 
a more rapid degradation in damping as the flutter point is approached. This 
is shown by the results given in figure 7 which compare the loci of mode 1 at 
M = 0.6, 0.7, 0.8, and 0.9. The tick marks represent calculations performed at 
dynamic pressure increments of 0.24 kPa. The results show that as Mach number 
is reduced the variations in the real part of the roots near the instability 
are increased for the same increment in dynamic pressure. Calculated flutter 
dynamic pressures and frequencies are given in table II. A value of equivalent 
viscous damping C of 0.005 was assumed for all calculations. 

Since modal damping is proportional to 

( Real part of root \ 

______ 

Imaginary part of root J 

the same analytical results presented as root locus plots can also be presented 
in the familiar form of damping and frequency versus dynamic pressure. An exam- 
ple of these data at M = 0.90 and M = 0.60 for modes 1 and 2 is presented 
in figures 8(a) and 8(b), respectively. These plots are useful in qualitatively 
assessing the nature of flutter onset. At M = 0.90 (fig. 8(a)) the damping of 
mode 1 is low throughout the dynamic-pressure range and the slope of the curve 
at flutter is not severe. These results indicate that the response in mode 1 
would be quite evident and the variation in frequency could be easily monitored. 
Figure 8(b) shows that at M = 0.60 the reduction in damping above a dynamic 
pressure of 6.7 kPa is quite rapid, which indicates a more violent flutter 
onset. These qualitative results were confirmed during tunnel tests. 
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System-off experimental flutter points were measured at M= 0.6, 0.8, 
and 0.9. Figure 9 presents a comparison of the predicted and measured flutter 
dynamic pressures and frequencies. Measured dynamic pressures and frequencies 
are given in table II. The results show good agreement at all Mach numbers. 


Wing With Active Flutter-Suppression System 

The eigenvalues of equation (9) are the roots of the characteristic equa- 
tion for the wing with active flutter suppression (system on) . Root locus 
plots at M = 0.90 and M = 0.60 are presented in figures 10 and 11. Root 
loci for each flexible mode and for the two filter modes introduced by the 
denominator terms (see eq. (10)) (s + 21.6)2 and (s 2 + 43.3s + 2922) are 
given in each figure. (It should be noted that the equations were solved 
with the complete transfer function defined by eq. (10). However, only the 
results of the 10 flexible modes and the 2 filter modes are presented in 
figs. 10 and 11.) At M = 0.90 an instability at a dynamic pressure of 
9.289 kPa results from a coupling between the first structural mode and the 
filter mode (s + 21. 6) 2 . An 84.8-percent increase in flutter dynamic pressure 
over that of the wing without flutter suppression is predicted. Calculations 
performed at M = 0.7 and 0.8 (not shown) predict the same type of behavior 
with increases in flutter dynamic pressure of 68 percent and 60 percent, 
respectively. Calculations at M = 0.60 (fig. 11) show a different behavior 
with the first flexible mode going unstable at a dynamic pressure of 9.337 kPa. 
This probably occurs because as Mach number decreases the system-off flutter 
frequency increases, which results in a decoupling between the first wing mode 
and the filter mode. Calculated system-on instability dynamic pressures and 
frequencies are given in table II. All system-on calculations include schedul- 
ing of the filter parameter Dj,, as given by equation (11). 

System-on tests above the system-off wing flutter boundary were performed 
only at M = 0.90. Figure 12 presents a summary of the predicted and measured 
effect of the flutter-suppression system. The model was stable to approximately 
42 percent above the system-off flutter boundary. At this point, the control- 
system-commanded aileron displacement exceeded that available on the model 
(6 a max = ±14°) and the system saturated, which resulted in an instability. 

The frequency of the instability was approximately 8.5 Hz, which indicated a 
complete loss in effectiveness of the flutter-suppression system when the sat- 
uration occurred. System-on tests at M = 0.8 and M = 0.6 were performed 
below the system-off wing flutter boundary. In all cases the wing was stable 
and the response of the wing to tunnel turbulence was reduced with the flutter- 
suppression system operating. 

It is believed that the saturation and resulting instability experienced 
on the model was a result of tunnel turbulence. Unpublished data taken from 
measurements of pressure fluctuations in the transonic dynamics tunnel indicate 
that the largest pressure peaks occur in the 8 to 1 5 Hz frequency range at Mach 
numbers between 0.87 and 0.95. At M = 0.90 and q = 7.59 kPa, equations (10) 
and (11) predict a control- surface displacement of approximately 11° per g unit 
for motion in the first flexible mode. At this test point accelerations in the 
frequency range of the first flexible mode were in excess of 1.6g. It can 
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therefore be assumed that the turbulence input to the model is greatest in the 
frequency range of the first mode and this led to large control-surface motions 
and finally to saturation. Prior to the instability, control-surface commands 
in excess of 20° were recorded. Reference 6 points out the need to consider 
stability and turbulence criteria simultaneously when designing a flutter- 
suppression system. 

In an effort to predict these results analytically, a gust calculation was 
performed on the wind-tunnel model using the methods described in appendix B. 
However, the power spectral density of the wind-tunnel turbulence is not modeled 
properly using a Von Karman gust spectrum. During the calculations the char- 
acteristic length of the Von Karman spectrum was varied to match the root- 
mean-square (rms) deflection of the control surface at the saturation point 
(M = 0.9; q = 7.59 kPa) . The measured and predicted variation in rms control 
deflection is given in figure 13. A gust intensity of 0.3048 m/sec was assumed 
for all calculations. The effect of reducing the gust length is to increase 
the gust input power at the first mode frequency. 

Since the design of a flutter-suppression system depends on both stability 
and response to turbulence, some effort to adequately model the tunnel turbu- 
lence is necessary if meaningful results are to be obtained. This study is 
beyond the scope of the present investigation. 


CONCLUSIONS 

An analytical method for predicting the increase in stability provided by 
an active flutter-suppression system has been presented. The method is based 
on approximating the unsteady aerodynamic forces in the time plane through an 
interpolating function in the frequency plane. The analytical method is 
applied to an aeroelastic model equipped with an active flutter-suppression 
system that was tested in the Langley transonic dynamics tunnel. Some of the 
important conclusions are: 

(1) The analytical technique presented provides a convenient method for 
adding active control systems to the equations of motion. 

(2) The use of interpolating functions to approximate the unsteady aero- 
dynamics provides a good prediction of the flutter characteristics of the wing 
without flutter suppression at all Mach numbers investigated. 

(3) Analytical results predict an 84.8-percent increase in the flutter 
dynamic pressure for the wing with flutter suppression at M = 0.90. Experi- 
mental results demonstrate a 42-percent increase prior to control-system 
saturation. 
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(4) Results of this study indicate the need for a more adequate model of 
wind-tunnel turbulence before a thorough evaluation of the analytical techniques 
can be performed for the wing with active flutter suppression. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
December 22, 1978 
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APPENDIX A 


AERODYNAMIC APPROXIMATION 

Three-dimensional unsteady forces are normally computed at a given Mach 
number for simple harmonic motion at specific values of reduced frequency k. 
The transfer function which relates control- surface motion to wing response is 
normally expressed as a ratio of polynomials in the variable s. This appendix 
describes a technique, similar to that described in reference 3, which permits 
the variation of the aerodynamic forces with reduced frequency to be approxi- 
mated by a rational polynomial in s. 

Consider the function 


Q(k) = A 0 + A-, (jk) + A 2 (jk) 2 + 


6_ A m (jk) 


m=3 (3k + 3^2) 

to be an approximate fit to Q. The real and imaginary parts of Q are 


(AT) 


Or = A o - A 2 k 2 + 




k 2 A, 


k 2 Ar 


k 2 A* 


k 2 + 3-j 2 k 2 + e 2 2 k 2 + b 3 2 k 2 + b 4 2 

B 4 kA 6 




(A2) 


B-] kA 2 B 2 kA 4 B^kAij 

Q x = A-| k + + + 


k 2 + B-| 2 k 2 + B 2 2 


k 2 + Bo 2 k 2 + B, 2 


J 


Q is calculated at discrete values of reduced frequency k. At each value of 
reduced frequency real and imaginary error functions are determined from equa- 
tions (A2) ; that is 


%,i = CR,i + lB R,i J * C > 
E I,i = Ql,i + lB Ifi J{c} 


where 


^R fi J = 




-v- 2 


-1 0 ki 


0 -ki 0 


.2 


-ki 2 -ki 2 


-k _• 2 


ki 2 + B] 2 ki 2 + B 2 2 ki 2 + B3 2 ki 2 + B 4 2 J 


— B t k i ^2^i ^ 3^i ~^4^i 


ki 2 + B, 2 ki 2 + B 2 2 ki 2 + B3 2 k £ 2 + e 4 2j 


(A3) 
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APPENDIX A 


{C} 


A 0 

A 1 






and i refers to a particular reduced frequency k^ at which Q is calcu- 
lated. Defining a complex error function as 

E i = E R,i + 3 E i , i 

a least-squares fit can be passed through the N data points by setting 


a N 

r 51 (Ei x E i*) = o 

3c i=1 

Where Ej* is the complex conjugate of E^. Performing this differentiation 
results in the following set of normal equations 


]EI ( Qr, i + lB Rfi J{C}){B Rfi } + (^Q lfi + lB Ifi J{C}){B Ifi } = 0 


which can be solved for the coefficients of the fit. That is, 

-|-1 

N / \ N v 

{C} = -\2~ UB Rfi }LB Rfi J + {B Z/i }lB Ifi J J QR,i{B Rfi } + Qi,i{B Ifi }j 

i=1 ' ' i=1 v ' 


(A4) 


Since k = iob/V, let s = jco and jk = s(b/V). Substituting this relationship 
into equation (A1 ) results in 


Q(s) 


Aq + A 1 


s + a 2 | 


b 

V 




A m s 



(A5) 


where the coefficients Ag, A-j, . . . , Ag are determined from equation (A4) . 
The values of 3 m _2 are arbitrarily selected from the range of reduced fre- 
quencies for which Q has been calculated. 

As described in reference 3, the form of equation (A5) permits an approx- 
imation of the time delays inherent in unsteady aerodynamics subject to the 
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APPENDIX A 


following requirements: complex conjugate symmetry, denominator roots in the 

left-half plane, and a good approximation of the complex aerodynamic forces at 
s = ju). The form of equation (A5) is used to fit all of the wing motion, con- 
trol surface, and gust unsteady aerodynamic forces. 
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APPENDIX B 


FORMULATION OF EQUATIONS OF MOTION FOR STABILITY AND GUST ANALYSIS 

The equations of motion are formulated in terms of real matrices by using 
an "approximation function" for the complex aerodynamic forces. The variation 
with s of the aerodynamic matrix Cq] is given as (see appendix A) 


[q] = Ia 0 1 + [A-, if - ]s + Ca 2 l( - ) s 2 + ]> 


6 [A m ls 


m=3 


r + l M 


After substituting [q] for 2 ([q]J in equation (4) , the equation of motion 
may be written as ' ' 


[Mis 2 + teSMWjJs + -PV 2 [q] + [K ] ) {q} = {f g } 
v 2 / 


(Bl) 


These are the equations for n structural modes with r active controls where 
Cm! represents the mass matrix, [r] the stiffness matrix, P the fluid density, 
V the fluid velocity, {Fq} the aerodynamic gust force, and {q} the response 
vector. All the matrices in equation (Bl) are of the size n * (n + r) . 


Flutter Analysis - No Controls 

Substituting the aerodynamic approximating coefficients into equation (Bl ) 
with {f g > = 0, the equations of motion in terms of real matrices are written as 


3 2[m] + [2£MW_]s + [r] {q> + -pv4 

J 2 




+ s(-Va,] + Ia 0 ] + 2“ 

W m=3 


[A m ] s 


s + ~ $m-2 

D 


{q} = 0 


(B2) 


vrtiere {q) = {q} for the no-control case. The matrices in equation (B2) are 
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APPENDIX B 


of the size n x n. Multiplying through by the denominator term yields a 
polynomial in s of the form 


(s)£[a] s 2 + [ 5 ] s + [cl] + A-J (s) [63 + A 2 (s)[e] + A 3 (s)[f] + A 4 (s)[G]j{q} = 0 


where 


(B 3 ) 


A 0 (S) = I s + l M( s + l Ml* + l + l ^ 


a-, (s) = sfs + £ e 2 j( s + ^3) f s + -■ P4J 


A 2 (S) = sfs + ^ 3J fs + £ 6 3 j Is + £ 3 4 j 


a 3 (s) = sfs + - 3 ) f s + - 3 2 ) fs + - 3 4 ] 


a 4 (s) = sf s + ^ 3J fs + ^ 3 2 j ls + l 3 3 j 


[a] = [m] + -pb 2 [A?] 
2 Z 


[b] = -pvb[A n ] + [2CMdJ 

2 


[Cl = — PV 2 [An] + [k 1 


[51 = -pv 2 [a 3 ] 
2 * 


[e] = -pv 2 [a 4 ] 
2 
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[p] 


-PV 2 [A 5 ] 
2 


[g! = ^pv 2 [a 6 1 


After the indicated polynomial products are performed, equations (B3) can be 
written as 

f[F 6 ls6 + [f 5 1s 5 + . . . + [F 0 ]) {q} = 0 (B4) 

The matrix coefficients tF^l (i = 0, 1, . . .,6) are functions of dynamic pres- 
sure and velocity for a given Mach number. By using the relationships that 


r 


S 5{qf 

s 4 {q> 


{X} = < . > 


^{q} J 


equation (B4) can be reduced to the following 6n first-order equations: 


s(x} = [a] {x} 

where 


[a! 


[-F 6 - 1 F 5 ] [-F 6 - 1 F 4 ] 

[i] 0 

0 [I] 


[-Fg-TF,,] 

0 

0 


0 


0 [i] 


0 


(B5) 


(B6) 
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The matrix [a] is 6n x 6n. For a fixed value of Mach number, dynamic pres- 
sure, and velocity, the eigenvalues of equation (B6) are the roots of the char- 
acteristic flutter equation. Since the matrix elements [f^] are functions of 
dynamic pressure (for a constant Mach number) , the loci of roots as a function 
of dynamic pressure can be constructed. 1 These loci correspond to the varia- 
tion in the eigenvalues of each flexible mode as dynamic pressure is varied at 
a constant Mach number. 


Stability Analysis - With Controls 

For the case of r controls, the response vector {q} can be expressed in 
terms of n structural modes and r control deflections as 





Equation (Bl), with {Fq} 


0, can be written as 



(B7) 


where the subscript c denotes a control quantity. The control law relates 
control- surf ace motion to wing response and can be written in the form 


{q c } = [T]{z t } (B8) 

where {z^-l are the values of wing response at the sensor location. The response 
{zj-} can be written in terms of modal response by 


{z^} = 


where [<|) t ] is the matrix of modal deflections at the sensor locations. 
Therefore, 


{q c } = [TH4> t Hq} 


TFor a given Mach number, the flight velocity V varies somewhat because 
of the change in speed of sound with altitude. For a wind-tunnel model, Mach 
number fixes the value of V. 
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Typically, the transfer function matrix [t ] is expressed as a rotational 
polynomial in s. Therefore, let 


[Tl 


D(S) 


(B9) 


where D(s) is a polynomial representing the common denominator of all the [t] 
terms, and [t n 1 is a matrix of the resulting numerators. (For the example in 
the text, t n = s 2 N(s).) Substituting equations (B9) and (B8) into (B7) results 
in the following equation of motion: 


[s 2 tM] + st2?Ma) n ] + [Kl]{q} + -PV 2 


/ b\ /b\ 

s 2 - [A 2 ] + s[ - ) [Ai ] + [A 0 ] 

\V/ \V/ 


6 [Ajjj] S 


1 

^ l \ 

{q} + 

s 2 [M c ] + [K c ] + -PV 2 

s + l ^ 


2 

\ b L 


L 


45 

\ V j 


A 2 ,q c 


-5' 

\VJ 


»q c 


A o,q c 


6_ [ A m,q c ] s 


m=3 


8 + b &m ~ 2 / 


[T N ] [<t. t ]{q} 


D (s) 


= 0 


where the values of ^A^ f q c j (i = 0, 1, . . .,6) are the aerodynamic matrix 

coefficients for each control surface. Multiplying through by the denominator 
term yields a polynomial in s of the form 


D(s)^A 0 (s)[[a]s 2 + [Bis + [C]] + A-, (s) [ 5 ] + A 2 (s) [e] + A 3 (s) [£] + A 4 (s) [G]) 

+ ^A 0 (s)[[A c ]s 2 + [B c ]s + [C c ] ] + A-, (s) [ 5 C ] + A 2 (s)[e c ] + A 3 (s)[F c ] + A 4 (s) [G c ]^ 


x [Tfl] [<|> t ] 


Kq} = 0 


(BIO) 
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where 


[A c l = W[A 2fq J + [M c l 
[B c ] = ^Pb[A 1/qc ] 


tc c ] 

[ 5 C 1 

ti c ] 

iFcl 

[G c ] 


Jpv 2 [Ao >qc ] + [K c l 

r 2 [ ft 5, qo ] 


Equation (BIO) can be written as 

^[F m ]s m + [Fm -1 ] s 111-1 + . . . + [Fq]^ {q} = 0 

where m = 6 + highest order of polynomial D(s). In a manner similar to that 
discussed in the previous section this equation can be reduced to a series of 
first-order equations of the following form: 


where 


s{ X> = [a c Hx} 


fs 


1,1-1 {q}] 


{X} 


s m-2{q} 

< • 


sP{q} 


(BIT) 
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and 



F m-1 






[aJ = 


[i] 

0 


0 

[I] 


(B12) 


0 0 [I] 0 


The matrix [ A c ] is mn x mn. For a fixed value of Mach number, dynamic pres- 
sure, and velocity, the eigenvalues of equation (B12) are the roots of the char- 
acteristic equation. Root loci can now be constructed which correspond to the 
variation in the eigenvalues of the system as dynamic pressure is varied. 


Gust Analysis - With Controls 

The gust response analysis is performed using power-spectral-density (PSD) 
techniques similar to those described in reference 7. Equation (B1 ) permits 
the direct evaluation of the system response to a sinusoidally varying gust. 
The term {^q} in equation (B1 ) is defined as 


{ f g } =--pv 


{ A 0 ,g} 




{ A m,g} s 



w 


g 


The modal response of the system with controls per unit gust velocity can be 
determined by solving the following set of simultaneous equations at discrete 
values of s (where s = jco ) . In this context the aerodynamic approximation 
is used only to interpolate data between calculated values of reduced frequency 


r • , 

[Equation (BIO)] = D(s) 

Wrr 


(<* 0 . 


A 0 < s ) I { A 0,g} + { A 1,g} s + 


a 2 , g } S 2 ] 


+ A-,(S){A 3f g} + A 2 (S) {A4 ^ g} + A 3 (S){A 5f g} + A4 ( S ) { Ag f g} 


(B13) 
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where 


{A 0 ,g} =-^PV{A 0 f g> 

{A 1/g } =-lpb{A 1/g } 

1 b 2 

{A 2,g> =--P -tA 2 f g> 

{A 3 , g } =-lpv{A 3fg } 

{A4 f g} =- ~P v ^ A 4,g-^ 

{ ^5,g> = ”PV{A 5( g} 

{A 6,g> ="PV{A 6#g } 

and {A^ f g} i = 0, 1, . . 6 are the aerodynamic matrix coefficients for a 

sinusoidal gust. The control-surface transfer function can then be evaluated 
by 


^ t^ {q} 

Wg D (S) Wg 


(B14) 


The power-spectral-density (PSD) values of control-surface motion are determined 
by evaluating 


4> 0 (U) = <fJ T (£D) | H((D) | 2 

where 

H(to) = Control-surface frequency-response function described by equation (B14) 
(j) j (oo) = Von Karman PSD gust spectrum defined by 
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L 


c 


+ 


8(7 .339I/JJ/V) 




,-iii/e 

+ (7 .339UJ/V) 2 J 


The root-mean- square (rms) value of control-surface motion per unit rms gust 
velocity cr w ^ is defined by 



(B75) 
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TABLE I.- FREQUENCY, GENERALIZED MASS, AND MODAL DEFLECTION DATA 


Mode 

Natural frequency, 
Hz 

Generalized mass, 

kg 

Modal deflection at 
sensor location 

1 

5.233 

3.678 

0.9228 

2 

19.129 

7.769 

-.6361 

3 

20.906 

7.044 

-.0002 

4 

25.769 

2.970 

.3450 

5 

46.110 

4.714 

.1760 

6 

61 .234 

4.758 

.2356 

7 

79.682 

5.156 

.0199 

8 

86.030 

11.297 

.0002 

9 

98.087 

7.558 

.0438 

10 

118.150 

5.501 

.0172 


TABLE II.- SUMMARY OF RESULTS 



Wing without flutter suppression 

Wing with flutter suppression 

Mach 

number 

Analysis 

Exper iment 

Analysis 

Experiment 

q» 

Frequency, 

q* 

Frequency, 

q» 

Frequency, 

q r 

Frequency, 


kPa 

Hz 

kPa 

Hz 

kPa 

Hz 

kPa 

Hz 

0.9 

5.027 

8.1 

5.36 

CO 

• 

o 

9.289 

4.1 

(a) 

— 

.8 

6.033 

8.9 

6.08 

8.6 

10.128 

4.2 

— 

— 

.7 

6.799 

9.6 

— 

— 

10.896 

4.0 

— 

— 

.6 

7.374 

10.3 

7.25 

10.1 

9.337 

7.6 

— 

— 


a No flutter to q = 7.590 kPa. 
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Figure 6.- Dynamic-pressure root locus at M = 0.90 (system off). 
Arrows indicate increasing dynamic pressure. 
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Figure 7.- Root loci of mode 1 as a function of Mach number and dynamic pressure (system off). 
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J Unstable 

■ 
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Imaginary part, rad/sec 



= 0.90 (system on) 
pressure. 
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Mode 6 
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Figure 11.- Dynamic-pressure root locus at M = 0.60 (system on). 
Arrows indicate increasing dynamic pressure. 
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Figure 12.- Effect of control law on flutter dynamic pressure 
as a function of Mach number. 
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